Reexamination of the long-range Potts model: A multicanonical approach 
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We investigate the critical behavior of the one-dimensional g-state Potts model with long-range 
(LR) interaction 1 /r d+a , using a multicanonical algorithm. The recursion scheme initially proposed 
by Berg is improved so as to make it suitable for a large class of LR models with unequally spaced 
energy levels. The choice of an efficient predictor and a reliable convergence criterion is discussed. 
We obtain transition temperatures in the first-order regime which are in far better agreement with 
mean-field predictions than in previous Monte Carlo studies. By relying on the location of spinodal 
points and resorting to scaling arguments, we determine the threshold value o- c (q) separating the 
first- and second-order regimes to two-digit precision within the range 3 < q < 9. We offer convincing 
numerical evidence supporting er c (q) < 1.0 for all q, by virtue of an unusual finite-size effect, namely, 
finite-size scaling predicts a continuous transition in the thermodynamic limit, despite the first- 
order nature of the transition at finite size. A qualitative account in terms of correlation lengths is 
provided. Finally, we find the crossover between the LR and short-range regimes to occur inside a 
narrow window 1.0 < a < 1.2, thus lending strong support to Sak's scenario. 
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I. INTRODUCTION 



Microscopic models with long-range (LR) interactions 
decaying as a power law, i.e., as l/r d+a , have aroused 
renewed interest during the last decade. Beyond their 
fundamental relevance to the understanding of critical 
phenomena, they have started playing a seminal role 
in the modeling of neural networks p] and spin glasses 
with Ruderman-Kittel-Kasuya-Yosida (RKKY) interac- 
tions |2j, systems undergoing phase separation, e.g., 
highly ionic systems Q and model alloys 0, and more 
widely in a large class of chemical or biological models 
where electrostatic interactions, polarization, or van der 
Waals forces play a central role. They have also attracted 
much attention in the framework of nonextensive ther- 
modynamics, where a possible equivalence with short- 
ranged (SR) models is under consideration 

Since the very early work of Ruelle [|| , LR spin mod- 
els in particular have been extensively studied. In one- 
dimensional models, it has been widely shown that long- 
range order occurs at finite temperature if and only if 
a < 1 0, Q, H, 0, , and this is in strong contrast to the 
SR case where no phase transition exists at finite tem- 
perature. Fisher and co-workers have shown that the 
upper critical dimension is reduced to d* — 2a, whereby 
one-dimensional LR models exhibit mean-field-like be- 
havior for a < 0.5, with the critical exponents taking 
on their classical values v = l/cr and 7 = 1 provided 
the phase transition is continuous. Conversely, the criti- 
cal behavior for a > 0.5 yields nontrivial exponents, and 
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LR models in effect go through a variety of universality 
classes as a is varied within this range, thus exhibiting 
rich critical behavior. Due to the ability to continuously 
vary the range of interaction, which in effect alters the ef- 
fective dimension of the model, one-dimensional LR mod- 
els are therefore a powerful paradigm for studying the 
dependence of critical properties on dimensionality, e.g., 
in systems above their upper critical dimension |l2j |. 

While significant emphasis has been placed on the 
Ising chain (see, e.g., [13| for a review), specific stud- 
ies of the LR g-state Potts model are less numerous and 
rather recent. These include a transfer matrix study com- 
bined with finite-range scaling (FRS) |l4j |. renormaliza- 
tion group (RG) analyses based on Wilson's momentum- 
shell method 0, or a real-space procedure 0, 0] , a 
cluster mean- field approach J19L and Monte Carlo (MC) 
simulations [H HJ ill HI UM, E3 The last, however, 
mostly focused on the case q = 3, and led to numerical 
estimates of critical exponents and temperatures showing 
some discrepancies. Due to the higher ground state de- 
generacy, this model reveals a phase diagram markedly 
richer than that of the Ising chain. It has been shown 
in the SR case that the transition turns from a contin- 
uous to a first-order one as the number of states q is 
increased beyond a threshold value q c (d) depending on 
the dimension ality of the model. For instance, q c {2) = 4 
and Q c (4) = 2 HaE3 ( see also Ref - E3 for a complete 
review) . As for the LR case in d = 1 , Glumac and Uzelac 
have shown from MC studies of the three- and five-state 
Potts model [2(j that the same sort of behavior occurs, 
i.e., there is a so-called tricritical point at some value 
<J c {q) depending on q, and the transition is continuous 
for a > a c . This qualitative picture was later reinforced 
in [23 for q — 3,5,7,9 and i n [ 241 for q = 3, both re- 
lying on MC studies, and in [21| using a graph- weights 
approach. On the other hand, it is noteworthy that RG 
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analyses dedicated to LR models have remained thus far 
rather inconclusive, where distinguishing between first- 
and second-order transitions is concerned [IsllTft ITU Il8l | . 

Although it is now believed that q c depends continu- 
ously on the range of interaction for this class of models, 
the exact location of the tricritical line separating both 
regions is still fairly controversial. The biggest hurdle 
for a precise and reliable determination of this border- 
line actually stems from the weakening of the discontin- 
uous transition as a c is approached, even if exceptionally 
large lattice sizes are simulated, e.g ., using the efficient 
Luijten-Blote cluster algorithm [23. While for q = 3 a c 
was claimed to lie between 0.6 and 0.7 in [2(j, Krech and 
Luijten pointed out that a — 0.7 still belongs to the first- 
order regime, and that the second-order regime may set 
in for (j = 0.75 only |24|. The situation with q = 5 turns 
out to be even worse, with numerical estimates available 
only within fairly large ranges: a lower boundary value 
of 0.8 was reported in [20], whereas 0.7 < cr c (5) < 1.0 
according to [22]. These results have not yielded a very 
precise phase diagram as yet, with the only reliable as- 
sertion being that tr c (g) increases with q. 

The marginal case a = 1 raises another set of thorny 
questions: in |22| it was reported that the phase transi- 
tion changes from a second-order to a first-order one for 
q > 9, while it has been shown b y K osterlitz, using a 
model with a continuum of states [30j, and later on by 
Cardy, using a discrete model ^t|> that inverse square 
interactions give rise to a Kosterlitz-Thouless (KTjtran- 
sition, i.e., one governed by topological defects [31|. It 
is worth mentioning that both hypotheses may be rec- 
onciled, at least partially, by following a scenario similar 
to the one devised in 

[H El IP, whereby for XY-like 
models with nonlinear nearest-neighbor interactions, the 
KT-like transition is preempted by a first-order transition 
whenever the nonlinearity becomes strong enough. While 
the recent work of Luijten and Messingfeld on the three- 
state Potts model [25j lends further support to Cardy's 
assertion, the controversy still appears unsettled, how- 
ever, and in this view a determination of the asymptotic 
behavior of cr c (q) as q — > 00 seems of major interest in- 
deed. 

We wish to shed light on some of these contradictory 
results using MC simulations in generalized ensembles, 
with particular emphasis put on the first-order regime. 
The aim of this work is thus twofold. First, we pro- 
pose an implementation of the multicanonical algorithm 
dedicated to the numerical study of LR models. This 
algorithm, devised by Berg and Neuhaus a decade ago 
[35L l36| . has been successfully applied in the past to SR 
models undergoing first-order transitions. As numeri- 
cal studies of models exhibiting first-order transitions are 
dramatically hampered by huge tunneling times when us- 
ing standard Metropolis update mechanisms [33, , a 
multicanonical approach is indeed an appropriate choice 
for both the determination of the location of the tricrit- 
ical line and the estimation of critical couplings in the 
first-order region a < a c . Our purpose is therefore to 



adapt the scheme initially proposed for SR models so as 
to make it suitable for a large class of LR models. Second, 
by relying on an extensive study for 3 < q < 9 and a wide 
range of a values, we arrive at convincing conclusions re- 
garding the location of the tricritical line, the range of 
validity of the mean-field-like behavior, which we find 
much larger than in previous studies, and the crossover 
from the LR to the SR regime, although the last was in- 
vestigated for the three-state model only. We show that 
our multicanonical implementation yields numerical esti- 
mates which are in agreement with and often better than 
those found in previous studies, although our simula- 
tions were performed by relying on medium lattice sizes, 
i.e., L < 400 spins. In particular, we obtain the follow- 
ing estimates for a c (q): er(3) = 0.72(l),cr(5) = 0.88(2), 
cr(7) = 0.94(2) and cr(9) = 0.965(20), and these results 
are highly precise. We also offer convincing evidence that 
the phase transition in the limiting case a = 1.0 is not 
of the first order for all values of q, by virtue of an un- 
usual finite-size effect. A detailed finite-size scaling (FSS) 
analysis conducted for q — 9 shows that, while the transi- 
tion belongs to the first-order regime at finite lattice size, 
its first-order nature wanes quickly enough as size is in- 
creased so that the transition tends to a continuous one in 
the thermodynamic limit. We give a qualitative account 
of this behavior in terms of correlation lengths, and by 
raising some open questions regarding the dynamics of 
first-order transitions in the LR case, we try to challenge 
the usual picture inherited from SR models. Finally, by 
relying on the shape of the specific heat and computing 
several moments of the magnetization, we conclude that 
a crossover between LR and SR regimes occurs inside a 
narrow window 1.0 < a < 1.2 



The layout of this article is as follows. In Sec. [HJ 
we first review some prominent features of the LR Potts 
model through a mean-field (MF) analysis. Special em- 
phasis is given to the calculation of the location of spin- 
odal points, a feature we will use in Sec. lIVI for estimating 
cr c (q). Section 11111 is devoted to implementation details 
of the multicanonical algorithm specific to LR models. 
We discuss the iteration procedure used to obtain the 
best estimate for the density of states, the choice of an 
efficient predictor, and a reliable convergence criterion. 
Improvements over the original algorithm are made in 
order to work out the algorithm instability due to low 
energy levels being unequally spaced. Numerical results 
regarding both first- and second-order regimes are then 
presented in Sec. lIVI Since we do not know of any previ- 
ous implementation of a generalized ensemble algorithm 
in the case of LR spin models, we pay particular attention 
to comparison with other standard MC algorithms, i.e., 
in terms of dynamical exponents, tunneling times, and 
accuracy of numerical estimates of critical couplings. 
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II. MODEL AND MEAN-FIELD THEORY 

Throughout this work we consider a ferromagnetic 
Potts model incorporating LR interactions in d — 1. This 
model is derived from a generalized q-state Potts Hamil- 
tonian, i.e., 



X J &(7i ,<Jj 



hidr, 



where the Potts spin Cj at site i can take on the values 
1, . . . , q, the first sum runs over all pairs of sites, and hi is 
an external aligning field favoring condensation in state 
(To- Incorporation of LR interactions is carried out by 
setting 

1 

Jij = J(\i-j\) = 



3 



\d+a ■ 



where d — 1 throughout this study, and a is an adjustable 
parameter which can be related to the effective dimen- 
sion of the model. As a falls off to —1, this model tends 
to the mean-field case where all interactions have equal 
strength, whereas the limiting case a — > oo corresponds 
to a pure SR model. Crossover from LR to SR behavior 
should actually take place at a — 1.0 [l^, 0], yet no 
numerical evidence has been given so far for this model 
which would reinforce this assertion. The thermodynam- 
ics of the model is studied numerically by way of the 
following order parameter: 



M 



q max„ p n 



where n = 1, . . . , q, and p n is the density of Potts spins 
in state n, which varies between 1/q at infinite tem- 
perature and 1 in the ground state. On a lattice of 
size L, numerical implementation is carried out by us- 
ing periodic boundary conditions, i.e., one adds up in- 
teractions between all the spins of the original lattice 
only, and replaces the bare coupling constant J(r) by 
J(r) = J2n=-oo J( r + nL). Retaining only interactions 
such that \i — j\ < L/2 leads indeed to strong shifts in 
energy and critical couplings for low a values, especially 
when finite-size scaling is to be used with medium lattice 
sizes. For the purpose of numerical evaluation, this sum 
may be reexpressed as 



J(r) = 



1 



1 



r l+cr 



1 



(7,1 



r 
L 



where £(s, a) denotes the Hurwitz zeta function. The 
self-energy will be omitted since it is just an additive 
constant to the total energy. 

Mean-field behavior can be readily obtained by using a 
variational MF approach (see, for instance, |41j). which 
relies on the minimization of the following functional 

F[p] = Tr pH + kT Tr p In p 



with respect to a trial density matrix p. Here the 
trace operation means a sum over all spin configurations, 
and the dependence of H and p on the spin configu- 
ration is implied. F[p] reaches a minimum whenever 
p = e~ H / kT /Z, i.e., in the case of a canonical Gibbs 
distribution, and this minimum yields the free energy 
of the system. The mean-field approximation allows us 
to express the density matrix p of the whole system as 
a product of one-site density matrices pi which depend 
solely on the spin variable at site i. We may thus rewrite 
the trace operation as a sum involving traces on single 
spin variables, namely, 

F[p] = Tr » Pi m Pi ~ X! h% Tr * piSa ' > CT ° 

i i 

- - ^2 Tr, pi Iv, i>,J, ,S, T , T . 

For further comparison with numerical results, we are 
mainly interested in expressing the free energy as a func- 
tion of an order parameter which is as similar as pos- 
sible to the one defined above. This is carried out by 
parametrizing the trial density matrix pi in terms of the 
following order parameter field: 



qS a 



1 



where the average is weighted by the trial density matrix 
Pi. Seeing that all states but state cro are equivalent, the 
constraint Tr pi — 1 thus yields 



p t (mi,<Ji) 



Considering a uniform external field hi = h, we have 
rrii = m for all sites; hence the free energy per spin f(m) 
reduces to 



q-\ 



hm — C(l + (i)m 2 + kT{(l — m) ln(l — m) 
l + m(q- 1) 



■ln[l + m(g-l)]} 



(1) 



where we dropped terms which are constant in m so that 
/(0) = 0, and £(1 + a) is the Riemann zeta function. 
This function expands as I /a around a = 0; hence tran- 
sition temperatures are expected to vary as 1/a in the 
vicinity of the MF regime. Equilibrium values of the or- 
der parameter are located at minima of the free energy, 
and it can be seen that m = is a stable minimum for 
kT > 2£(1 + a)/q. For q = 2, there is no third-order 
term in the series expansion of f{m); hence a second- 
order transition occurs at kT c — £(1 + a). For q > 3, the 
negative coefficient in the third-order term of the series 
expansion creates a second minimum, which physically 
corresponds to a first-order transition. At the transi- 
tion temperature, the free energy has the same value at 
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FIG. 1: Reduced temperatures of spinodal points fcTi/£(l+<r) 
and kT-ijC,{\ + a) together with the reduced transition tem- 
perature fcT c /£(l + a), as a function of q in the MF approxi- 
mation. 

both minima. Following |28j |. the exact transition tem- 
perature kT c may be computed by simultaneously solving 
/(to) = /'(to) = and yields 

kT c _ q-2 
C(l + cr) ~ (q-1) In(q-l)- 

Similarly, spinodal points are computed by jointly solv- 
ing /'(to) = /"(to) = 0, giving temperature points at 
which either one of the two minima vanishes. These 
equations possess one trivial solution, namely, kT\ — 
2£(1 + cr)/q corresponding to the extrema at m — be- 
coming unstable, and a nontrivial solution kT2 which may 
be obtained numerically by solving the following equa- 
tion, 

2 q-1 \ V 2 y 

where S = 1 + yj\ + 2(1 - q)/{Kq) and we have set 
K = ((1 + a)/kT2- Alternatively, one may also express 
/ as a function of the MF energy E = — £(1 + <j)m 2 
and impose f'(E) = f"(E) — 0. While these equations 
yield the same kT\ and kT2 as above, the two expres- 
sions of / obviously do not have the same shape. Spin- 
odal points are sketched in Fig. ^ f° r Q between 2 and 
10. These correspond to the limit of metastability for 
each subphase, respectively. For temperature points ly- 
ing inside this temperature range, there exist two values 
of the order parameter corresponding to a null curvature 
of the free energy, a feature which is known to induce 
a long-ranged (i.e., low wave number) instability. This 
in turn triggers a phase transition through the so-called 
spinodal decomposition [42|. As expected, the width of 
the spinodal curve T2 — T\ shrinks to zero as q — > 2, and 
accounts for the second-order nature of the transition at 
q = 2, since in this limit the two minima merge into 
a single large minimum responsible for the well-known 
divergence of fluctuations at a continuous transition. 



III. THE MULTICANONICAL ALGORITHM 

The Metropolis algorithm (hereafter denoted as be- 
longing to the class of canonical algorithms, i.e., rely- 
ing on a Boltzmann weighting) has long been considered 
the paradigm for Monte Carlo simulations in statistical 
physics, yet this method faces some severe drawbacks in 
situations where the sequence of states created by the 
Markovian chain leads to very repetitive dynamics, i.e., 
dramatically low acceptance rates and exponentially di- 
verging autocorrelation times: this makes it necessary 
to simulate systems over exceedingly long runs in order 
to obtain good statistics and reliable estimates of ther- 
modynamical averages (see, for example, 1 13) and the 
contribution by Krauth in 44] for an introductory re- 
view). This is the case when one comes to simulating 
systems with rugged free energy landscapes, e.g., poly- 
mers, proteins, and disordered systems including spin- 
glasses, for the dynamics may then get trapped in one 
of numerous local minima, especially at low tempera- 
ture. One experiences similar behavior when simulating 
first-order phase transitions (the so-called supercritical 
slowing down H3), where the tunneling time between 
coexisting phases grows exponentially with the system 
size, due to the increasingly high free energy barrier to 
be overcome (e.g., 0). 

Since slow dynamics mainly results from the combina- 
tion of weighting the Markovian chain with Boltzmann 
weights and using local updates, there have been several 
attempts to devise efficient update algorithms based on 
global updates, e.g., cluster algorithms, which in the case 
of continuous transition decrease critical slowing down 
by several orders of magnitude (see 0, 0] ; also a LR 
implementation in | 2 9j |). On the contrary, multicanonical 
methods j3j| l37l l48j are based on random walks in 
the energy landscape, irrespective of the particular move 
update utilized, whereby aflat energy distribution is now 
sampled. First, this results in the algorithm quickly sam- 
pling a much wider phase space than in the canonical 
case, by allowing the system to cross any free energy 
barrier. Second, this allows the density of states to be 
computed over the whole energy axis, thus extending the 
reliability of reweighting procedures over a much wider 
range of temperature than in the case of standard his- 
togram methods |49j , where poor histogram sampling at 
low-energy usually induces str ong statistical bias. As op- 
posed to multihistogramming j5fj , a single run is needed 
to cover the energy range of interest. Once a reliable esti- 
mate of the density of states has been obtained, it is then 
straightforward to compute thermodynamical functions 
otherwise hardly within reach of canonical simulations, 
e.g., canonical entropy and free energy. It is notewor- 
thy that this simulation technique actually belongs to a 
larger class of algorithms called generalized- ensemble al- 
gorithms, which encompasses variants based on random 
walk in the entropic variable (" 1/fc ensemble" or "en- 
tropic sampling" algorithms |5lll52| '). or the temperature 
variable (e.g., "simulated tempering" |U|3|). 
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Rationale 



B. Iteration scheme 



The rationale behind the multicanonical algorithm is 
the generation of a Markovian chain of states {<Tj}, whose 
weights W mu {E(ai)) are tweaked so that one eventually 
gets a flat energy histogram, i.e. if P(E) denotes the 
probability in energy and n(E) is the density of states, 

P mu (E) oc n(E)W mu (E) = const. 

Since n(E) usually increases drastically with energy, low- 
energy states are thus sampled much more often than 
high-energy ones. 

Following Berg in [Hlf. we compute W mu (E) through 
an iterative procedure, starting from an initial canoni- 
cal simulation at inverse temperature /3q. /3q indirectly 
sets the energy below which the energy histogram is to 
be flat, i.e. E max = (E)^. Thus, fcTo = l//3o must 
be chosen high enough to ensure that the final energy 
histogram spans a suitably large energy range upward, 
e.g., reaches the energy of the disordered phase in the 
case of a first-order transition, and extends even further 
away if one wants to observe with satisfying accuracy the 
free energy plateaus signaling the limit of metastability. 
For convenience, we now define an effective Hamiltonian 
H m u(E), so that 

W mu (E,p ) = e-P° H ^ E \ 

Hence, multicanonical simulation can be envisioned as a 
canonical simulation at inverse temperature /3g with the 
usual Boltzmann weight, provided the original Hamilto- 
nian is replaced by an effective Hamiltonian to be deter- 
mined iteratively. As a side note, a cluster implementa- 
tion in the framework of the multicanonical algorithm is 
thus far less straightforward, since this effective Hamilto- 
nian has fundamentally a global nature, whereas canoni- 
cal simulations explicitly preserve the locality of the orig- 
inal Hamiltonian (see, e.g., the multibond approach in 

Si). 

Denoting H^ U (E) as the true estimate of the effective 
Hamiltonian, we may thus write 

n(E) oce^""^. 

The microcanonical inverse temperature /3(E) may be 
easily related to H^ U (E), as we have (assuming k = 1) 

m = ^l = ^ u (E) 



dE 



dE 



Since the dynamics of the Markovian chain is gov- 
erned by the transition rate W(a — > b) = 
min(l,exp{f3o[H mu (E a ) - H mu (E b )]}), we may write, 
for two states infinitely close in energy, i.e., whenever 
E b = E a + 5E, W(a ->• b) = mm(l,exp[-/3(E a )SE}). 
Hence it is the microcanonical temperature which is the 
relevant quantity where the dynamics (e.g., the accep- 
tance rate) of the multicanonical algorithm is concerned. 



We initially set H^ U (E) = E, or equivalently /3°(E) = 
(3q, as this indeed corresponds to a canonical simulation 
at temperature l//3o- At step i, a simulation is performed 
using a Boltzmann weight with effective Hamiltonian 
H l mu (E)\ then an energy histogram N l (E) is eventually 
computed using independent samples. Incidentally, tak- 
ing truly independent samples proves useful during the 
late stages of the iteration scheme only, where the aim is 
then to refine a nearly flat histogram. During early iter- 
ation steps, histograms may be computed using noninde- 
pendent samples without significantly affecting the con- 
vergence. We now denote E l min as the lowest energy level 
that was reached throughout the previous runs, includ- 
ing step i: this is the energy level below which H^(E) 
will have to be predicted, since no histogram data are 
available inside this energy range. Issues regarding ade- 
quate predictor choice will be considered later on in this 
section. The rules for updating H^} at step i + 1 from 
P % mu at step i are based on the following equations. For 
E > E max , H^(E) — E, i.e., the dynamics is canonical- 
like at inverse temperature (3q for all iteration steps. For 



P'rain — E < Emax, 



(3 l+1 (E) = (3 l (E) + 



where 



9b 



go ln N l (E- 
SE 



SE) 



N*(E) 



(2) 



and <7q is a raw inverse damping factor proportional to 
the reliability of the fcth histogram. It has been shown 
in |55j . following an error calculation argument, that 



9o 



N(E)N(E + 5E) 
N(E)+N(E + 5E) 



provides an estimator proportional to the inverse of the 
variance of (3 i+1 (E). Once (3 l+1 (E) is known, H l +^(E) 
is derived by a mere integration starting from the initial 
condition H mu (E max ) = E max . Finally, for E < E % min , 
Hmu (E) wm have to be computed using a suitably cho- 
sen predictor, until at last E l min becomes equal to the 
ground state energy. A cubic spline is then fitted to 
H mu (E) at every bin center, and this curve is used to 
compute acceptance probabilities during the next run. 

It can be seen that Eq. leads to a steady state when- 
ever N(E) is constant over the energy range of interest. 
Writing a recursion equation involving /3(E) instead of 
H mu (E), together with the inclusion of a damping factor, 
allows one to handle the situation where some bins have 
null entries, a case which otherwise leads to a fairly spiky 
graph for H mu (E) and inconsistent dynamics. Accidental 
null entries at energy values E or E+SE will simply leave 
(3(E) unchanged, and the corresponding parts of H mu (E) 
thus move as a block. Since acceptance rates hinge on 
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FIG. 2: Lowest energy levels for q = 5, a — 0.5, N = 400, 
computed by sorting energy samples from a long simulation 
run. Each level is drawn as a horizontal line. 



the microcanonical temperature, this in effect drastically 
reduces bias on the dynamics. Considering a small set 
of histogram bins that are copiously filled for the first 
time during a given iteration run (e.g., high-energy bins 
during the early iteration runs whenever we begin with a 
canonical simulation) , we see that the related cumulative 
inverse damping factor first soars and produces a great 
amount of change in (3(E) in the couple of runs that fol- 
low, and then decays progressively to zero as these bins 
continue to be filled. By taking into account all the data 
that have been sampled up to step i, this modified re- 
cursion both clearly stabilizes the algorithm and reduces 
relative errors due to poor histogram sampling. 

Choosing the most appropriate value of the histogram 
bin width results from a trade-off between resolution and 
computation time. A higher resolution on the one hand 
guarantees good histogram flatness, and is especially cru- 
cial at low energy levels, where the density of states dis- 
plays a rugged graph. On the other hand, we impose a 
fixed number of independent samples per histogram bin, 
so as to give the histogram variance an acceptably low 
value; hence a low SE implies more simulation steps per 
iteration. Our approach is thus to first choose a fairly 
high SE, e.g., one yielding around 20 bins, during the 
early stages of the iteration process in order to obtain 
a rough picture of the density of states, and then to 
progressively reduce SE once the ground state has been 
reached. As will become obvious in Sec. IIV Al the ulti- 
mate value of SE deeply affects the attainable precision 
on the computation of spinodal points, since the latter is 
based on a precise location of free energy plateaus, and 
this indeed entails having enough bins belonging to a 
given plateau. As a rule of thumb, the best compromise 
is then to obtain between 100 and 300 histogram bins 
in the final stage, with the number of bins increasing as 
the a value corresponding to the second-order regime is 
approached. 

In this view, the unequal spacing of energy levels in LR 
spin models deserves specific attention. As witnessed in 



Fig. [21 large energy gaps separate isolated energy levels 
or tiny groups thereof in the vicinity of the ground state, 
whereas the distribution gradually turns into a near con- 
tinuum above E ~ —1025. Setting a low SE value leads 
in turn to nonaccidental null entries in those bins located 
inside energy gaps, whereby (3(E) never gets updated at 
isolated energy levels and go is always zero. Since the 
graph of the density of states looks indeed fairly wrinkled 
near the ground state, and the dynamics there is notice- 
ably sensitive to even the smallest departure of H mu (E) 
from the ideal line, we would then observe a sharp steady 
peak in the lowest part of the energy histogram, which 
the present recursion would not be able to suppress. One 
could trivially think of working this out by implementing 
variable- width bins that would span energy gaps. This 
is, however, impracticable since the distribution of energy 
levels is not known prior to starting the iteration process 
(for this is precisely what we intend to compute with 
the density of states). To circumvent this limitation, we 
have modified the previous recursion so that null entries 
are always skipped, however accidental or nonaccidental 
they may be. Denoting by E a and E bl with E a < E b , the 
centers of histogram bins located on each side of a set of 
contiguous empty bins, we have 



(3 i+1 (E a )=(3 i {E a ) + 



gi(E a ) ln NHE b ) 



E b - E a N i (E a ) 



(3) 



where (3(E a ) = (3 {H mu (E b ) - H mu (E a )} and we now 
impose 



go(Ea) 



N(E a )N(E b ) 
N(E a ) + N(E b )' 



hence go can never be zero. In order to avoid losing 
details of the shape of H mu (E) for E a < E < E b that 
were possibly collected during previous runs, we update 
H mu (E) through a linear difference scheme, 

SH mu (E) = 5H ^ Eb J ~ 5 H mu (E a ) {E _ Ea)+SHmu{Ea): 

&b — Ea 

where 5H mu (E) = H*£(E) - H l mu (E). While this has 
obviously no effect where nonaccidental null entries are 
concerned, this favors quicker convergence during the 
early runs where the inadequate shape of H mu (E) is more 
likely to produce empty bins. 

The iteration process stops whenever the energy his- 
togram has become suitably flat over the energy range 
of interest, namely, between the ground state energy and 
E m ax for our purpose. We evaluate this property by com- 
puting the standard deviation of histogram entries, as 
well as the same quantity for the logarithm of histogram 
entries restricted to nonempty bins. The latter seems to 
be a better indicator since it is sensitive to both poorly 
populated bins and histogram peaks, whereas the former 
increases only with rather spiky histograms. In addition, 
we estimate the degree of convergence of the algorithm 
by computing the mean square distance between HJ mu (E) 
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and H^(E) after the ground state has been reached. 
We then compute a threshold value for each indicator by 
trial and error, based on a couple of short runs for various 
lattice sizes and bin widths. 



C. Reweighting procedure 

Once H mu (E) has been satisfactorily computed, a long 
production run is performed using this effective Hamil- 
tonian in place of the original one, and then estimates 
of thermodynamical quantities of interest at inverse tem- 
perature (3 are computed using a reweighting scheme, i.e., 
formally, 



{A) = 



Y.e (A) E n{E)e 



-0E 



where (A) E denotes the microcanonical average of A 
at energy E, and the partition function is given by 
Z = J2e n(E)e~ l3E . The best estimate for the density 
of states n(E) is provided by n(E) oc N(E)e PoHmu(E \ 
where N(E) stands for the number of bin entries at en- 
ergy E computed from the production run. In order to 
avoid numerical overflows, as well as to suppress bias 
resulting from possibly strong variance on microcanoni- 
cal averages, we found it more appropriate to compute 
(A)g from a sum running over samples instead of en- 
ergy bins, i.e., (A) = Y Jl A i w ( E i)/Y J i w ( E i)i where 
w(Ei) = e PoHm U {Ei)-PEi-K _ K jg then determined so 
as to avoid both numerator and denominator overflows. 
Providing that the histogram sampled during the produc- 
tion run is flat to a good approximation, the maximum 
m e p o H mu (E)-0E ig reac hed whenever dH mu (E)/dE ~ 
/?//3o) which yields the energy value at which K is to be 
computed. In addition, since the reweighting scheme in- 
volves an exponential contribution of H mu (E), the result- 
ing curve e @oH mu (E)-t3E j g gangly peaked around the 
maximum; hence it is clear that only histogram points 
in the vicinity of this maximum contribute to (^4)^- In 
effect, we found that the existence of two distinct max- 
ima, or equivalcntly of two energy values for which /3(E) 
has the same value, coincides with the occurrence of a 
first-order phase transition. 

Following the same reweighting procedure, we compute 
partial free energy functions, i.e., F(j3,m) where m is 
the order parameter, and reweighted histograms of the 
energy, i.e., N rw (j3,E). The partial partition function is 
straightforwardly derived from a partial sum over sam- 
ples having the prescribed order parameter, 



(4) 



which then yields F(j3,m) — — \nZ(j3,m)/f3. Similarly, 
a reweighted histogram of the energy is obtained from 



D. Predictor choice 

We now discuss some issues related to the choice of an 
efficient predictor for E < E m i n . For small lattice sizes, 
we initially feed the algorithm with an effective Hamil- 
tonian H mu (E) = E, and the objective is then to find 
an appropriate trade-off between speeding up the con- 
vergence of E l min toward the ground state and avoiding 
algorithm instability. While the former demands that 
H mu (E) have a sufficiently high slope below E min , the 
latter still requires that the algorithm remain ergodic to 
a suitable extent. Our implementation relies on a first- 
order predictor, H mu (E) = a + bE, and we impose conti- 
nuity on H mu (E) at E m i n . The simplest approach is then 
to choose a predictor slope so that continuity on H' mu (E) 
is enforced at E = E min , i.e., b = (3(E min )/ (3 . While 
E m in reaches the ground state rather quickly using this 
predictor, the dynamics often gets locked in very low en- 
ergy levels due to the particularly steep slope of H mu (E) 
in the vicinity of the ground state. The time needed by 
the iteration scheme to recover from this deadlock and 
obtain a flat histogram thus becomes prohibitive. On 
the other hand, choosing 6=1 leads to the smoothest 
yet slowest convergence, and avoids deadlock issues. An 
efficient compromise is thus to ensure a "weak" continu- 
ity at Emin, i-e., by computing the slope of the predictor 
using a least-square scheme based e.g., on the first 10% 
of points above E m i n . 

For large lattice sizes where reaching the ground state 
energy can become time consuming, we resort to a "scal- 
ing trick" whereby H mu (E) is initially guessed from the 
density of states obtained at a smaller lattice size. This 
approach was initially mentioned by Berg and Neuhaus 
|37| , and claimed to work perfectly within the framework 
of a study of the two-dimensional ten-state Potts model 
with nearest-neighbor interactions where the energy is 
additive to a perfect extent. The presence of LR interac- 
tions, however, slightly worsens the case, especially at low 
a. The scaled density of states is computed as follows. 
Let us consider, for the sake of simplicity, two systems E 
and E with respective lattice sizes L = N and L = 2N, 
and let us divide the latter into two subsystems Ei and 
E2 of equal size L. Since H mu (E) — kTohin(E), where 
n(E) stands for the density of states, we have to compute 
n(E) for S as a function of n(E) for E. Neglecting the 
interaction between subsystems Ei and E2, and denoting 
by Ei the energy of Si, the density of states for E just 
reads n(E) ~ J2 El n(Ei)n(E — Et), which yields 



)] 



5E J 



where SE is the energy histogram bin width. Provid- 
ing that n(E) is a monotonic and rapidly increasing 
function of E, we may use a saddle-point approxima- 
tion to evaluate the former sum. The maximum of 
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FIG. 3: Dots indicate the initial guesses Hmu(E) that were 
fed into the iteration scheme at L = 400, q = 5, and a = 
0.3(0), 0.5(+), 0.9(D). Each initial guess was computed using 
Eq. i.e., by scaling a true estimate obtained at L = 200. 
Solid lines show true estimates H mu (E) as obtained after the 
whole iteration scheme at L = 400 converged. The straight 
dashed line sketches the original Hamiltonian, i.e., H mu (E) = 
E. 
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FIG. 4: Energy histogram as computed after indicated runs, 
for q — 5, a = 0.9, L = 400 spins, using Eq. © to compute 
the initial effective Hamiltonian H mu (E) from a previous run 
at L = 200 spins. Labeling on y axis indicates normalized 
probabilities. 



H mu {E\) + H mu (E — E\) is reached whenever E\ = E/2; 
hence 



H rnu (E) ~ 2H mu I — I + kT In — (5) 

This expression may be readily extended to lattice sizes 
that are any multiple of the original size. Figure |3| 
sketches results obtained for q = 5 and a = 0.3,0.5, and 
0.9. A series of iteration runs was first conducted with 
L = 200 spins in order to obtain an estimate of H mu (E) 
for this lattice size, then this estimate was scaled using 
Eq. (JSJ and used as the initial guess H m u (E) for the next 
series of iteration runs at L = 400. Equation (J5J yields a 
very acceptable guess for a = 0.9, and the two curves are 
hardly distinguishable from each other. As illustrated in 



Fig. the energy histogram becomes nearly flat within 
five iterations. For a = 0.3 and 0.5, the agreement re- 
mains quite satisfying, yet the initial guess falls slightly 
below the true estimate at low energy levels, and the 
lowest-energy bins are exceedingly enhanced during the 
first iteration runs. More iteration runs are thus required 
to obtain a perfectly flat histogram as a is decreased, and 
the benefit of this approach in effect becomes negligible 
for er < 0.3. Indeed, the algorithm then spends a great 
number of iteration steps being trapped in low energy lev- 
els, seeking to rectify the shape of the density of states in 
this energy region until convergence is obtained: starting 
from an initial canonical effective Hamiltonian actually 
yields better performances. Since, for systems with LR 
interactions, computation time scales with L 2 , using this 
"scaling trick" thus greatly reduces the time needed for 
proper convergence, at least for a > 0.3, and partially 
compensates for the lack of a hybrid multicanonical- 
clustcr algorithm dedicated to LR models. 



E. Algorithm performance 

In order to measure the performance of our implemen- 
tation, we have computed a dynamical exponent z de- 
fined as the scaling exponent of a relevant characteristic 
time t of the simulation, i.e., t oc L z , where L denotes the 
lattice size: while for second-order transitions it is widely 
known that the integrated autocorrelation time repre- 
sents such a relevant time, for first-order transitions the 
tunneling time through the energy barrier (rtun) proves 
to be a more meaningful indicator [3£|. We define the 
latter as one half of the average number of Monte Carlo 
steps per spin (MCS) needed to travel from one peak 
of the reweighted energy histogram (N rw (f3, E)) to the 
other, with the temperature being set to the transition 
temperature. Tunneling time is expected to grow expo- 
nentially with L for canonical algorithms, and to scale as 
a power law of L for multicanonical algorithms Is?! . In 
both cases, it appears that the chosen characteristic time 
is a good indicator of how quickly the demands in CPU 
time should grow with increasing lattice size: for second- 
order transitions, this is the time needed to generate 
truly independent samples, while for first-order transi- 
tions, this tells us at what rate the dynamics spreads out 
over the energy barrier and thus to what extent samples 
get efficiently picked from the two phases in coexistence. 

The integrated autocorrelation time is computed by 
using the well-known time-displaced correlation function 
which displays an exponential-like short-time behavior, 
namely, < I ) mm (t) ~ e~*/ T ; r is then derived from a simple 
integration scheme. Since the latter function makes sense 
within equilibrium only, we first discard n thermalization 
steps, where n is obtained by using the nonlinear relax- 
ation function that describes the approach to equilibrium 
|43j and averaging over several dry runs. An interesting 
point regarding multicanonical simulations is that, since 
they are random walks in energy space, "thermalization" 
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FIG. 5: Integrated autocorrelation time r vs lattice size L for 
q = 7 and a — 0.2, 0.4, 0.6, 0.8. Dynamic exponents computed 
from a fit to L z are z = 1.09(1), 1.15(1), 1.38(1), 1.55(1), re- 
spectively. 
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FIG. 6: Tunneling time T% un vs lattice size L for q — 7 and 
a — 0.2, 0.4, 0.6, 0.8. Dynamic exponents computed from a fit 
to L z are z = 1.25(1), 1.30(2), 1.37(1), 1.53(1), respectively. 



(although this term is no longer appropriate as far as 
generalized ensemble algorithms are concerned) always 
occurs rather rapidly; simulations based on a nearly flat 
histogram have shown that a value of 1000 MCS is indeed 
appropriate on average. 

Results for q — 7 and a lying between 0.2 and 0.8 are 
shown in Fig. [S] for integrated autocorrelation times, and 
in Fig. for tunneling times. The slight dispersion in 
the power-law fits arises from the fact that simulations 
at larger sizes were conducted with a higher number of 
MCS between measurements in order to reduce mem- 
ory overhead. Where computing tunneling times is con- 
cerned, this results in some tunneling events being pos- 
sibly skipped and the average tunneling time being over- 
estimated. Both figures show, however, that a power-law 
behavior is perfectly plausible. In the case of first-order 
transitions, the reduction in simulation costs is thus dras- 
tic in comparison with standard canonical algorithms. 

For both indicators, we obtain an average z slightly 
above 1.0 for a = 0.2, yet z increases smoothly with 
decreasing range of interaction. This may be accounted 
for by the fact that spatial and time correlations grow 
as we depart from the MF regime and approach the SR 
one. As for tunneling times, the prefactor turns out to 
be slightly higher near the MF regime, and z increases 
at a lower rate with increasing a than is the case for 
autocorrelation times. 

Since there are no other numerical studies of LR mod- 
els based on multicanonical algorithms to our knowledge, 
comparison is limited to estimates obtained for SR mod- 
els. For the three-state Potts model, canonical simu- 
lations using local updates led to z — 2.7 [5(|, while 
Swendsen and Wang obtained z ~ 0.6 using their per- 
colation cluster algorithm ^(|. For further comparison, 
the Metropolis algorithm applied to a SR Ising chain in 
d = 2 and d = 3 yielded a value of z slightly above 2 
[53, whereas Wolff's cluster algorithm led to z ~ 0.27 
|58| . While our value is slightly greater than in the case 



of cluster implementations, it is worth underlining that 
our multicanonical implementation yields reliable statis- 
tics within a single MC sweep, whereas several are needed 
in the case of a standard canonical simulation, whatever 
reweighting procedure may be used. 

IV. NUMERICAL RESULTS 

We have conducted multicanonical simulations for q G 
[3,9], using for each value of q an appropriate set of a pa- 
rameters between 0.3 and 0.9, so that we could observe 
strong and weak first-order transitions, as well as contin- 
uous ones. For q = 3, some simulations were performed 
with a > 1.0 in order to investigate the crossover from 
LR to SR regimes. Once the density of states had been 
determined using the iteration process described above, 
a production run was performed for lattice sizes between 
L = 50 and L = 400. The number of MC sweeps needed 
for each production run was computed so as to yield ap- 
proximately 5 x 10 4 truly independent samples. In this 
view, rapidly increasing autocorrelation times in effect 
restricted our work to lattice sizes L < 400. 



A. Free energy functions and FSS 

As already stated in the Introduction, a precise de- 
termination of the tricritical value a c {q) is a real chal- 
lenge, due to the weakening of the first-order transition 
as <7 C is approached from below. This makes traditional 
indicators e.g., latent heats or energy jumps, fairly inef- 
ficient, since observing clear jumps in the vicinity of the 
tricritical value entails simulating huge lattices. Glumac 
and Uzelac in 20] used three less traditional estima- 
tors, namely, the interface free energy, the specific heat, 
and the reduced fourth-order Binder cumulant, which all 
turned out to be less sensitive to this weakening: in par- 
ticular, the last quantity defined as Ul = (E 4 ) / (E 2 ~) 
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is expected to reach a nontrivial constant Uoo ^ 1 as 
L — > oo at a first-order transition only |5£j; by extrapo- 
lating to the thermodynamic limit from measures taken 
at different lattice sizes, they found a c to fall between 
0.6 and 0.7 for the three-state model. Still and all, this 
approach imposes simulating fairly large lattices (around 
L = 3000) for the extrapolation procedure to be reliable, 
let alone the fact that Binder cumulants may experience 
uncontrollable crossover effects |60j. Due to the modest 
lattice sizes that are within reach of our local update al- 
gorithm, we rather resort to an approach based on the 
location of spinodal points, which may be accurately de- 
termined already for medium lattice sizes. In marked 
contrast to multihistogram techniques, the multicanoni- 
cal method indeed allows one to obtain partial free energy 
functions (or, equally, reweighted histograms of the en- 
ergy) over a range of temperature which extends fairly far 
away from the transition temperature, with remarkably 
modest numerical effort. 

The basis of our method relies on the fact that the 
temperature difference between the two spinodal points 
will tend to zero as a c is approached, since there are no 
metastable states in the case of continuous transitions. 
Stated differently, the conditions under which metasta- 
bility occurs, i.e., both first and second derivatives of the 
partial free energy are zero, are met only at the critical 
point for a continuous transition: hence metastable states 
merge into a single large minimum as the first-order tran- 
sition turns into a second-order one. Such behavior has 
indeed been widely observed, e.g., for liquid- vapor tran- 
sitions near the critical point, and is borne out by our 
MF calculation. 

For a given set of (q, a) parameters, we determine 
the location of the spinodal points by first comput- 
ing a partial free energy function of the order param- 
eter [F(kT, m), see Eq. 0] over a large temperature 
range. Alternatively, we make use of a similar function 
of the energy, i.e., F e (kT,E) = - lnN rw (kT, E), where 
N rw (kT, E) denotes the reweighted histogram of the en- 
ergy. While the latter function plays the same role as 
the partial free energy of the magnetization, it yields a 
higher precision at low q, as we will witness in a moment. 
The limit of metastability at finite lattice size is then 
defined by dF e /dE = d 2 F e /dE 2 = 0, or alternatively 
dF/dm = d 2 F/dm 2 — 0: for a first-order transition, this 
condition is met at two temperatures T\ and T2 which 
satisfy the inequality T\ < T c < T2, where T c denotes 
the transition temperature. 

Since these free energy functions usually have rather 
rugged graphs, we first filter out rapid oscillations by 
means of a linear smoothing filter, whose order is com- 
puted so that we are left with at most three extrema over 
the whole temperature range of interest. By continuously 
varying kT within this range, we determine the temper- 
ature of each metastable state by monitoring the change 
in the number of minima (see Fig. UJ. In contrast to 
|20| . the transition temperature T C (L) is then obtained 
by imposing that the number of bin entries in N rw (E) 
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FIG. 7: Graphs of F e (kT, E) = - In N rw (kT, E) for q = 
5, a — 0.3, N = 400, at four characteristic temperatures: 
T\,T2,T eqh , and T eqw — T c denote the temperatures of the 
two metastable states, and the temperature of equal peaks 
heights, and that of equal peak weights, respectively. E/L 
denotes the energy per spin. 
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FIG. 8: Average energy per spin for q — 3, a — 0.2, L = 400, 
computed over both phases (solid line), ordered phase only 
(lower dashed line), and disordered phase only (upper dashed 
line). Vertical dotted lines indicate the four characteristic 
temperatures: from left to right, lower limit of metastabil- 
ity (fcTi), transition temperatures (equal heights kT eq h, then 
equal weights kT eqm ), and upper limit of metastability (kT^). 



be the same below and above the energy corresponding 
to the maximum of F e (kT,E). This corresponds to the 
so-called equal- weights condition as proposed by Lee and 
Kosterlitz in |60j , and is equivalent to the condition that 
the average energy be the arithmetic mean of the en- 
ergy of each phase. For the sake of comparison, how- 
ever, we also compute the temperature T eq h(L) at which 
both minima of F e (kT, E) have the same value. We then 
proceed with the computation of similar quantities us- 
ing F(kT,m), and we estimate statistical errors using a 
bootstrap procedure. 

Graphs of the free energy F e (kT,E) in Fig. [3 show 
that the peak and the plateau corresponding to the disor- 
dered phase are much narrower than those of the ordered 
phase. As a result, the precision in the determination 
of the temperature T-y of the lowest metastable state is 



11 



fairly lower than that of the upper metastable state (I2). 
This asymmetry increases with increasing q, and in effect 
precludes the use of reweighted histograms for the esti- 
mation of spinodal points at q > 7. For q = 9, we thus 
relied on the extrema of the partial free energy F(kT, m), 
since this function then becomes nearly symmetric and 
displays peaks that are well separated. Incidentally, the 
asymmetric shape of F e (kT, E) can be accounted for by 
the fact that specific heats have a different magnitude 
in the disordered and ordered phases, since this thermo- 
dynamic quantity is simply proportional to the standard 
deviation of the associated Gaussian peak [5!| . This may 
be readily observed by reweighting thermodynamical av- 
erages over a single phase at a time, once the maximum 
of F e (kT,E) which separates the two phases has been 
located. Figure |S] shows how this procedure was applied 
to the computation of the mean energy per spin of each 
subphase for q = 3, a — 0.2, and L = 400 spins. A simple 
visual inspection allows one to assess a much lower spe- 
cific heat for the disordered phase than for the ordered 
phase. 

At finite lattice size, however, all these temperatures 
experience a distinct shift proportional to the distance 
from the thermodynamic limit. Assuming that the FSS 
theory developed in for first-order transitions is also 
valid in the LR case, we therefore compute tempera- 
tures at infinite lattice size by assuming power-law cor- 
rections in 1/L. We also expect temperatures defining 
the limit of metastability to obey the same scaling be- 
havior, although the phenomenological theory proposed 
in |59j does not explicitly handle them. The inclusion 
of a second-order term proves necessary in order to to 
obtain satisfying fits, due to the presence of small lat- 
tice sizes in our set of data. Yet, interestingly enough, 
fitting finite-size temperatures to a power law of the 
form T(L) = T(oo) + aL b yields very similar extrapo- 
lated values, with discrepancies smaller than 0.1%, i.e., 
within our range of uncertainty. In addition, we observed 
that F e (kT,E) and F(kT,m) led to distinct finite-size 
shifts, with the latter function easily allowing one to drop 
second-order correction terms without much affecting the 
final result. 



B. Transition temperatures 

For the sake of completeness, we also compute tran- 
sition temperatures by relying on two other estimators, 
namely, the magnetic susceptibility, which for magnetic 
systems has more pronounced peaks than the specific 
heat, and Binder cumulants of the magnetization defined 
as C/W = 1 - (to 4 ) /(3(m 2 ) 2 ). The latter are known 

to cross at a critical fixed point defining the true 
critical temperature, yet, since the crossing point drifts 
smoothly over our range of lattice sizes, we assume a 
power law of the form L w for U^(L), with an unknown 
exponent w |6l| . In addition, these two quantities are 
advantageously used to obtain critical temperatures in 
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FIG. 9: Spinodal curve for 0.3 < a < 0.8 (q = 5). The tran- 
sition temperature T c is indicated by filled squares, and the 
limits of metastability T\ and T2 by triangles and diamonds, 
respectively. Errors are smaller than the size of symbols, and 
lines are drawn to guide the eyes. The dotted line shows the 
transition temperature as predicted by MF theory. 
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FIG. 10: Partial free energy F(kT, m) for q = 9, a = 0.3, L = 
400 (solid line) , together with the MF prediction (dashed line) 
as given by Eq. Q. 



the second-order regime as well (see Sec. IIV El for more 
details on this issue). 

Results for all temperature estimates are summarized 
in Table for q = 3,5,7,9, and sketched in Fig. for 
q = 5. As expected according to FSS theory, both defini- 
tions of the transition temperature, i.e., using equal peak 
weights vs equal peak heights, lead within error bars to 
the same estimates at infinite lattice size. Other quanti- 
ties T c ( x ) and T C (U^) yield very similar results, with a 
discrepancy never exceeding 1%. 

For all values of q, the transition temperatures pro- 
gressively depart from the MF line as a is increased. For 
5 = 5, for instance, the ratio between T c (x) and the MF 
value ranges from 97.3% at a = 0.3 to 83.9% at a = 0.8. 
We further notice that, for a given range of interaction, 
the adequacy of MF results is clearly improved at high q. 
As illustrated in Fig.UHlfor q = 9, a = 0.3, and L = 400, 
this agreement also holds, even at finite lattice sizes, for 
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TABLE I: Estimates of the critical temperature in the first- and second-order regimes (the latter is indicated by an asterisk): MF, 
mean-field; x> using location of peaks of the susceptibility; t/' 4 ' using crossing points of Binder cumulants of the magnetization; 
eqw,eqh, using the free energy, where T c corresponds to equal peak weights and heights, respectively; Ref. [20|. MC study based 
on multihistogramming and the Luijten-Blote cluster algorithm (g = 3) and a standard metropolis algorithm (q — 5); Ref. 
[T^ l. cluster mean-field method combined with an extrapolation technique based on the VBS (Vanden Broeck and Schwartz) 
algorithm; Ref. 13), transfer matrix method combined with FRS. 



q 


a 


T c (MF) 


Tc(x) 


T C (U (4) ) 


T c (eqh) 


T c (eqw) 


T c (Ref. [20]) 


T c (Ref. [19]) 


T c (Ref. [14]) 


3 


0.2 


4.034 


3.97(1) 


3.98(1) 


3.94(1) 


3.97(1) 


3.70° 




3.7023 




0.3 


2.836 


2.72(1) 


2.72(1) 


2.71(1) 


2.71(1) 


2.70 a 


2.71669 


2.5893 




0.4 


2.240 


2.086(4) 


2.089(6) 


2.075(5) 
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1.31638 
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0.5 
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1.052(2) 


1.051(2) 
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0.7 
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0.793(2) 
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0.794(2) 
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0.705(2) 


0.704(1) 


0.704(1) 
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"Refers to l/K e (AF). 
'Refers to l/K E {U^). 



the shape of the partial free energy F(kT, m) and the 
position of metastability plateaus. For q — 3 and q = 5, 
we can readily compare the transition temperatures with 
earlier MC studies. Results obtained in [2(j using either 
the Luijten-Blote cluster algorithm (q = 3) or a standard 
metropolis algorithm (q — 5) are in fairly good agree- 
ment with ours within an error bar that does not exceed 
1%, except in the case a — 0.2, where our estimate lies 
much closer to the MF prediction. We further compared 
our estimates with those obtained in |19| using a cluster 
mean-field method, and in [Tlj using a transfer matrix ap- 
proach. As illustrated in Tabled results obtained using 
the cluster mean-field approach combined with the VBS 
extrapolation algorithm yield a perfect match, with a de- 
viation as low as 0.1% on average over the whole range 
of er values. The discrepancy with estimates obtained 
using the transfer matrix method is slightly higher and 
amounts to 2% on average, except for low values of a 
where the agreement of our results with the MF predic- 
tion is, here again, far better. 



C. Change of regime 

As can be viewed in Fig. [§| spinodal points merge 
slightly above a ~ 0.8 for q = 5, and this indeed sig- 
nals a change of the nature of the transition. By plotting 
dkT m = kT2—kTi against 1 /a, we observe that for all val- 
ues of q the points fit quite well on a line for low enough 
a, and the slope of this line tends toward that of the MF 
curve. The case q = 7 is sketched in Fig. ^] where it is 
clear that the point at a = 0.6 marks the border between 
the linear and nonlinear behavior, illustrating the weak- 
ening of the first-order transition as a c is approached. 
Since temperatures appear to scale as 1 j a in the vicinity 
of the MF regime, it is thus more appropriate to work 
with T\/T c and T2/T c , for the scaling factors will then 
cancel out neatly except when approaching a c (q). As 
mentioned above, the latter ratio, which is sketched in 
Fig. 1121 offers a higher precision through a larger free 
energy plateau. As a falls off to the MF regime, this 
ratio tends, within error bars, to the value predicted by 
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FIG. 11: Difference between temperatures of metastability 
dkT m = kT2 — kTi vs 1/a for q = 7 (circles connected by 
solid lines). Errors are smaller than the size of symbols. MF 
prediction is shown for comparison (dashed line). 



FIG. 13: Phase diagram computed using FSS properties of 
spinodal points, for a < 1.0. Dotted lines are shown to guide 
the eyes. 
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FIG. 12: T 2 /T c vs 1/a for g = 3,5,7,9. Solid lines indicate 
polynomial fits. Dotted lines are guides to the eyes. Error 
bars are smaller than the size of dots, except where explicitely 
indicated. 

the MF theory, i.e., T 2 /T c = 1.01, 1.037, 1.059, 1.077 for 
q = 3, 5, 7, 9, respectively. On the leftmost side of the 
graph, we witness a sharp decrease of TijT c as a — > a c . 
This brings a quite reliable way of determining cr c (q) 
without much ambiguity, as opposed to, e.g., methods 
using the interfacial free energy or Binder cumulants. 
By fitting data points to a polynomial of degree 2 for 
q = 5,7, 9, and of degree 3 for q — 3, which turned out 
to yield the lowest error, we obtained the following nu- 
merical estimates: 

q c c 
3 0.72(1) 
5 0.88(2) 
7 0.94(2) 
9 0.965(20) 

The graph of a c (q) is sketched in Fig.^|for convenience. 
Considering the global shape of this graph, it is reason- 
able to expect a c (q) — > 1 as q — > oo. This would be 



clearly consistent with Cardy's scenario (as mentioned 
in the Introduction), according to which the border case 
a = 1.0 corresponds to a KT-like transition governed by 
topological defects. 



D. Unexpected FSS behavior of correlation lengths 
and the dynamics of first-order transitions 

Let us now briefly inspect the case q = 9, a = 1.0, 
where a simple analysis based on the shape of the free 
energy at a given lattice size might be markedly mislead- 
ing. In |22f , a first-order transition for q > 9 was reported 
on the basis of the observation of a double-peaked energy 
histogram. We have performed a series of simulations at 
L = 50, 100, 150, 200, 300, and 400 for this set of param- 
eters and computed corresponding (finite-size) spinodal 
temperatures Ti(L) and Tb(L) using the partial free en- 
ergy F{kT, m). As may be noticed in Fig. O a striking 
feature of this limiting case is the existence of metastable 
states at all finite lattice sizes, with a first-order character 
strongly enhanced at low sizes, despite the fact that FSS 
theory yields T 2 — T\ = in the thermodynamic limit. 
It turns out that the transition is clearly not of the first 
order in the thermodynamic limit, and this feature was 
also confirmed for q = 6,7, and 8; for q < 6, a precise 
location of metastable states became impracticable. 

At first blush this behavior significantly contradicts the 
expected picture, whereby for first-order transitions, the 
correlation length is finite and roughly independent of 
the lattice size, and is merely connected to the size of 
clusters. As a result, such transitions appear as if they 
were continuous until the lattice size overtakes the corre- 
lation length. With regard to SR models, this has been 
the standard scenario thus far, yet we feel strongly that 
this scenario may be somewhat challenged, at least qual- 
itatively to begin with, where models incorporating LR 
interactions are concerned. 

To set the stage for an attempt to interpret this be- 
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havior, we first turn to the consequences of finite lat- 
tice size on long-wavelength fluctuations when simulat- 
ing LR models with algebraically decaying interactions. 
The key point in the following discussion is the nature 
of the phase transition as observed from numerical data 
obtained at finite lattice size. On a lattice of size L 
with periodic boundary conditions, the largest allowed 
distance between any two spins is L/2, and this also cor- 
responds to the smallest interacting potential affordable 
on a given lattice. It is obvious that these spins experi- 
ence a stronger interacting potential whenever L is small, 
and hence the whole array of spins may be rigidly tied 
to an adequate extent for an order-disorder transition to 
occur through metastability. When increasing the lattice 
size, on the contrary, spins being a distance L/2 apart 
now experience weaker interaction, and this results in a 
softening of the transition. Whether this softening might 
be sufficient to yield a change of nature of the transi- 
tion at some (either finite or infinite) lattice size, so that 
the transition may be continuous in the thermodynamic 
limit, is however an unsettled question; this assumption 
is borne out at least for q = 9 and a = 1.0, as witnessed 
by our results. Alternatively, we may say that the trun- 
cation of LR interactions at small lattice size artificially 
shifts the model toward the MF regime, since the inter- 
acting potential now varies smoothly over the available 
distance of interaction. 

As seems obvious to us, the usual physical meaning 
attributed to the correlation length in the case of SR 
models, i.e., roughly speaking the average size of a clus- 
ter of contiguous spins having the same value, may no 
longer hold in the case of LR models: since all the spins 
of the lattice, however distant they may be, are tied to- 
gether through an interacting potential, there is basically 
no need of a long-range order for two distant spins to al- 
ready have slightly correlated fluctuations. In the context 
of first-order transitions, this means that either clusters 
may extend well beyond the size permitted by the value 
of the correlation length, or the correlation length itself 
may become infinite in the thermodynamic limit. This 
behavior has indeed already been reported in models of 
DNA thermal denaturation as well as in the context 
of wetting [6^| . 

In addition, we have performed simulations in the first- 
order regime at the finite-size transition temperature. We 
used, however, a Metropolis algorithm, since the associ- 
ated dynamics is closer to the real nucleation or spinodal 
decomposition picture than with a multicanonical algo- 
rithm. We observed indeed that clusters in the ordered 
phase always spanned the entire lattice, whatever the 
lattice size. As soon as the dynamics jumps from the 
disordered to the ordered phase, which we monitored by 
comparing the energy with the location of reweighted en- 
ergy histogram peaks, a single cluster forms very rapidly 
and nudges its way through the crowd of disordered spins 
so that it swiftly occupies the whole lattice. Thus, if both 
phases coexist insofar as, e.g., the energy histogram has 
a double-peaked structure, they actually do not coexist 
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FIG. 14: Linear fit of finite size temperatures vs 1/L for q = 
9, a — 1.0. Dotted, solid, and dashed lines correspond to fcTi, 
kT eq h, and kT2, respectively. Error bars are smaller than the 
size of symbols, except where explicitly indicated. In the limit 
L — > oo, the difference between temperatures of metastability 
tends to 0.0012. Within our error bars, the transition is thus 
clearly not of the first order. 



at the same time and merely alternate in time, as op- 
posed to what is considered the usual SR picture. In this 
respect, we would like to raise some challenging ques- 
tion regarding the dynamics of first-order transitions in 
the LR case: (i) Since both phases do not coexist at the 
same time, what physical meaning should be given to the 
interfacial free energy? (ii) Does a mechanism similar to 
nucleation take place in a LR system, and if in the affir- 
mative, how can it be reconciled with the mechanism of 
cluster growth involved in SR models? 



E. Beyond the tricritical line: From LR to SR 
behavior 



We now focus on some critical properties in the second- 
order regime o~ c (q) < a < 1.0, then we investigate the 
crossover from LR to SR behavior. As mentioned in |64| . 
"standard" FSS theory is valid for LR systems provided 
the effective upper critical dimension d* — 2a is greater 
than the geometrical dimension d = 1, i.e., a > 0.5. Thus 
for q > 3 we assume "standard" finite-size scaling equa- 
tions to be valid. We first determine the critical expo- 
nent v using nth-order cumulants of the magnetization, 
i.e., V n — d\n (m n ) / dj3, which have minima obeying the 
scaling law V™ m oc L l v [6{|. Our approach is to com- 
pute two numerical estimates of v by fitting reweighted 
averages of V™ m and V™ m to a power law of the lattice 
size, and then to average over both values. Other criti- 
cal exponents, i.e., /3 and 7, are computed using similar 
scaling laws, i.e., M(T c (oo)) oc L~ /5 / l/ , and X max « L~</ u . 
Figure El shows a power-law fit of peaks of V\ , Vi and 
X against the lattice size obtained for q = 5, cr = 0.9. 
Points lie neatly on a straight line when using a log-log 
scale, and give the following estimates: l/v\ = 0.668(2), 
l/i>2 = 0.669(2), j/i> = 0.940(4). Error bars were com- 
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scale, for L = 50, 100, 150, 200, 400 (q 
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TABLE II: Critical exponents in the second-order regime a > 
a do), and q = 3, 4, 5. Shown for comparison are results from 
Ref. U4[ obtained using a transfer matrix method, and from 
Ref. | 22| using a MC histogram approach. 
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(3/v 
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0.8 


0.624(6) 


0.574 


0.842(5) 


0.101(5) 




0.9 


0.54(1) 
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0.908(5) 
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1.0 






0.96(1) 


0.025(8) 


4 


0.8 


0.71(1) 
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0.940(4) 
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1.0 






0.97(1) 


0.04(1) 




1.0 [22] 






0.966 


0.017 



puted using a bootstrap procedure. Once v is known, we 
fit finite-size temperatures T C (L) defined from peaks of 
the magnetic susceptibility to a power law of the form 
T C (L) = T c (oo) + \L~ l / v and obtain an estimate of the 
critical temperature. With regard to critical couplings 
obtained from Binder cumulants of the magnetization, 
we follow the same procedure as in the first-order regime. 
Finally, the critical exponent f3 is determined by fitting 
A/(T c (oo)) to a power law of the lattice size, and slowly 
varying the temperature at which M is to be sampled 
until the best fit is obtained. In the example considered 
above, this leads to (3/v — 0.103(2). Results for other 
pairs of (q, a) values are summarized in Table [I] and Ta- 
ble [H] For the borderline case a = 1.0, only exponent 
ratios are shown. It can be seen that our estimates match 
fairly well those obtained from a previous MC study (2^ , 
and that the discrepancy with results obtained from a 
transfer matrix approach in 0] never exceeds 8%. As 
opposed to the conjecture made in |18|. the exponent v 
does clearly depend on q. 

If the relation a = 2 — r\ derived in |ll| is indeed exact 
for q > 3, we should thus observe the simple behavior 
7/1/ = 2 — r) = a in the second-order regime. As il- 



lustrated in the fifth column of Table [nj the qualitative 
behavior follows the conjecture, yet clearly a < 2 — 77, 
and the discrepancy is remarkably higher for q = 5 than 
for q = 3. Moreover, while it appears to shrink to as 
a — ► 1, it is unclear whether j/u varies linearly with a, 
considering the small number of points available. 

In order to get a deeper insight into the crossover to 
the SR regime, we then conducted several simulations 
at q = 3 for a above the borderline value a co = 1. This 
value has been reported to play the role of a critical range 
of interaction beyond which a crossover from LR to SR 
behavior sets in. According to 0,H3, o co — 2 — tisr, 
where t]sr denotes the value of the 77 exponent in the 
SR case. Since 7/^=1 for all values of q in the SR 
case, rjsR = 1, and this indeed leads to <j co = 1. It 
should be noted, however, that this definition, as initially 
proposed by Sak in [3^ on theoretical grounds, as well 
as the exact location of a co within the interval [1.0,2.0], 
is still controversial. As shown in Table [HJ 7/1/ indeed 
appears to reach its SR value as a — » l - , yet this ratio 
proves no longer reliable above the borderline value, as we 
will witness in a moment, and reliance on other quantities 
becomes necessary. 

We first review some exact results concerning the SR 
regime, which we obtained using an exact transfer matrix 
method. For q = 3, the transfer matrix is a 3 x 3 matrix 
having three eigenvalues, which in zero external field read 
Ai = 3cosh(/3/2) - sinh(/3/2), A 2 = A 3 = 2sinh(/3/2), 
where /3 = X/kT. By retaining the largest eigenvalue Ai 
only, and taking the limit L — > 00, we successively obtain 
the free energy per spin 



and the specific heat 

a 08) = 



ln(2 + e?) 

2/3 2 



(sinh/3/2-3cosh/3/2) z 



From there on, the correlation length is then computed 
using the standard formula £ = l/ln(Ai/A2), which then 
yields 



in 



3coth/3/2 - 1 



Finally, the magnetic susceptibility is obtained using the 
fluctuation-dissipation relation, which gives 

It is then straightforward to show that 
lini/3^00 lnx(/3)/ ln£(/3) = 7/^ = 1. However, eval- 
uating this ratio at finite inverse temperature, i.e., 
for a finite correlation length as imposed by a finite 
lattice size, yields a greatly overestimated result. For 
instance, we obtain 7/1/ ~ 1.3 for L = 400, a feature 
which is supported by our simulation results, e.g., 
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FIG. 16: Specific heat for various lattice sizes and a — 
1.0, 1.1, 1.2, 1.7 and the pure SR case (obtained using a trans- 
fer matrix method), from right to left. Data for other values 
of a have been omitted in order to preserve the clarity of the 
figure. C v was computed using the fluctuation-dissipation 
relation C v = ({E 2 } - (E) 2 )/ (kT 2 L). Solid, dashed, dot- 
ted and long-dashed styles refer to L = 50, 100, 200, and 
400 respectively, except for the SR case where they refer to 
L = 5, 10, 100, 200. 



7/V = 1.02(1), 1.14(1), and 1.23(1) for a = 1.1, 1.5, and 
4.0, respectively. Since the last two values are clearly 
overestimated, this in effect indicates the presence of 
exponential divergences and as a by-product drastically 
slow convergence of the correction to scaling. 

This analysis was corroborated by a study of the shape 
of the specific heat, which turns out to provide the most 
tractable approach at medium lattice sizes where distin- 
guishing between the SR and the LR regime is concerned. 
In the thermodynamic limit, C v {f3) admits a maximum 
C raax _ QJ618 at kT rn = 0.3767. It is enlightening to 
investigate the nonmonotonic behavior of this maximum 
at finite L, and this may be carried out by computing 
F(f3, L) and then C v {(3, L) while retaining all three eigen- 
values. Since the calculation is fairly involved, and the 
final result admits no simple expression, we shall here- 
after simply refer to the corresponding curve sketched 
in Fig. El When L is increased, the peak of the specific 
heat first increases to a maximum, and then graphs of C v 
collapse and merge gently as the thermodynamic limit is 
approached. Whenever it is witnessed in graphs obtained 
from simulation data, this feature thus signals a SR-like 
behavior. 

Simulations were performed for 1.0 < a < 4.0 for var- 
ious lattice sizes between L = 50 and L = 400, and we 
set the initial canonical temperature to fcTo = 1-0 so 
that the maximum of C v would be clearly visible within 
the whole range a > 1.0. As appears obvious from a 
glance at Fig. ED the cases a = 1.0 and a — 1.1, on 
the one hand, and a > 1.2, on the other hand, display 
fairly distinct qualitative behaviors. For a = 1.0, the 
specific heat reaches its maximum monotically, at least 
for the lattice sizes that were investigated. The slow- 
ing down in the increase rate as 1/L — > allows one 
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FIG. 17: Maximum of the specific heat vs inverse lattice size 
for a = 1.0, 1.1, 1.2, 1.3, 1.5, 1.7, 2.0, 3.0, 4.0 from top to bot- 
tom. The solid line is a reminder for the (exact) SR case in 
the thermodynamic limit. Other lines are guides to the eyes. 




FIG. 18: Magnetization vs kT for a = 1.0, 1.1, 1.2, and 1.7 
from right to left. Solid, dashed, dotted and long-dashed 
styles refer to L = 50, 100, 200, and 400 respectively. 



to assess a finite maximum in the thermodynamic limit, 
and this clearly shows that C v is a nondivergent quan- 
tity, thus bringing support to Cardy's scenario whereby 
the transition has a KT-like nature on the borderline 
a = 1.0. The same behavior is observed for a = 1.1. On 
the contrary, the qualitative behavior is clearly different 
for <7 > 1.2, where the maximum of C v first decreases 
with increasing lattice size, and then quickly reaches a 
plateau reminiscent of the exact SR behavior investigated 
above. While this plateau only slowly reaches the exact 
SR value as a — + 4.0 (see Fig. I17|) . we can however con- 
clude that the behavior is already SR-like. This assertion 
can be further confirmed by considering the magnetiza- 
tion, as sketched in Fig. El Graphs of this quantity 
clearly merge slightly above m = 0, whenever a > 1.2; 
hence there is no transition at finite temperature. While 
for a = 1.1 there remains some ambiguity due to statis- 
tical errors, for er = 1.0 the curves now clearly intersect 
around kT ~ 0.7, which at least shows that the behavior 
is no longer SR-like. We finally compute critical temper- 
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atures from the crossing points of Binder cumulants of 
the magnetization. We obtain j3 c = 3.3, 6.5, and 19 for 
(7 = 1.1, 1.3, and 1.5. As for a = 1.7 and a = 2.0, cumu- 
lants no longer cross except at kT — within statistical 
error (the latter case yielding /3 C between 150 and 200, 
yet with excessive dispersion). While the crossover ap- 
pears to take place in the very vicinity of the borderline 
a = 1.0, the critical temperature actually dies off quite 
slowly to as a increases. 

All these numerical results lend support to Sak's sce- 
nario for a > 1.0, namely, that a crossover from LR to SR 
behavior occurs whenever a co — 2 — t)sr. Nonetheless, it 
is worth mentioning that we found this crossover to occur 
within the finite, yet narrow range 1.0 < a < 1.2, and the 
pure SR case to be reached in the limit a — > oo only. We 
feel strongly that this is consistent with the RG scenario 
of Theumann and Gusmao [16j, whereby the crossover 
actually results from a competition between SR and LR 
fixed points. This competition, as seems obvious to us, 
may not resolve instantly whenever a crosses the border- 
line, and may thus blur this borderline over some finite 
region. 

V. CONCLUSION 

We have studied some critical properties of the long- 
ranged Potts model using a multicanonical implementa- 
tion of generalized ensemble algorithms. Our implemen- 
tation of the iteration procedure needed to obtain the 
density of states was shown to yield satisfying estimates 
of this quantity over a large range of energy and with 
much quicker and more stable convergence than with 
the initial historical algorithm. The multicanonical al- 
gorithm allows one to efficiently circumvent the slowing 
down traditionally experienced at first-order transitions, 
and at the same time makes the reweighting approach a 
fairly straightforward way of examining thermodynamic 
quantities over a large range of temperature with strik- 
ingly modest numerical effort, i.e., by simulating over 
medium lattice sizes and performing a single long simu- 
lation run. We have used this multicanonical approach 
to locate spinodal points in the first-order regime over a 
large range of q and a parameters. The shape of the spin- 
odal curve in the vicinity of the change of regime then 
yielded precise estimates of the tricritical value cr c (q) up 
to two digits. In particular, the value <r c (3) = 0.72(1) 
is perfectly consistent with the lower bound of 0.7 pro- 
posed by Krech and Luijten |24|. yet in terms of precision 



this is markedly better by an order of magnitude. In 
this respect, our multicanonical implementation allows 
us to obtain numerical results whose accuracy is at least 
comparable to that of previous numerical studies based 
on multihistogramming and the LR cluster algorithm, al- 
though our simulations were performed on lattices having 
fewer than 400 spins. We feel strongly that this approach 
might be successfully applied to other spin models incor- 
porating LR interactions, e.g., continuous spin models or 
frustrated systems. 

In addition, our study significantly extends the range 
of available estimates of critical couplings and exponents. 
In the first-order regime, the agreement with MF predic- 
tions, and in particular with Tsallis's conjecture T c ~ 1/er 
in the limit a — > |5J], is exceptionally good. In the 
second-order regime, the relation rj = 2 — a, conjectured 
to be exact for q — 2, is shown to yield an increasingly 
high discrepancy when q is increased, and its validity may 
just be reinforced in the vicinity of a — 1.0. We found 
however that the crossover from the LR to the SR regime 
occurs between a = 1.0 and a = 1.2, thus lending strong 
support to Sak's conjecture. Our detailed FSS analy- 
sis of the case q = 9, a = 1.0 yielded one of the most 
surprising results of this study, namely, the unexpected 
behavior of correlation lengths whereby the transition ap- 
pears to be of the first order at finite lattice size, despite 
the fact that FSS theory predicts a continuous transition 
in the thermodynamic limit. We feel strongly that this 
may be accounted for by the truncation of the LR po- 
tential, which artificially brings the model closer to the 
MF regime, yet we also pointed out that the physical 
meaning of the correlation length should be somewhat 
challenged in the case of LR models. The exact nature 
of the transition in the borderline case a = 1.0, however, 
needs further investigation, especially at large q where no 
results have been made available thus far. In this view, 
an efficient combination of a global update scheme with 
a multicanonical approach would be of prior importance 
to reach far higher lattice sizes. 
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